Introduction to Bayesian pharmacometric data analysis using NONMEM
Hands-on 4: PopPKPD analysis using MCMC
Metrum Research Group, LLC
August 19, 2019
Population Parameter Estimation
Estimation of Friberg-Karlsson model of drug-induced myelosuppression.
Specification of prior distributions
Setup R environment
rm(list = ls())
gc()
modelName <- "friberg3"
scriptName <- paste(modelName, "Rmd", sep = ".")
fitModel <- FALSE
## Relative paths assuming the working directory is the script directory
## containing this script
scriptDir <- getwd()
projectDir <- dirname(scriptDir)
dataDir <- file.path(projectDir, "data")
modelDir <- file.path(projectDir, "model")
outDir <- file.path(modelDir, modelName)
figDir <- file.path(projectDir, "deliv", "figure", modelName)
tabDir <- file.path(projectDir, "deliv", "table", modelName)
invisible(dir.create(figDir, recursive = TRUE))## Warning in dir.create(figDir, recursive = TRUE): '/data/IntroBayesNONMEM/
## deliv/figure/friberg3' already existsinvisible(dir.create(tabDir, recursive = TRUE))## Warning in dir.create(tabDir, recursive = TRUE): '/data/IntroBayesNONMEM/
## deliv/table/friberg3' already exists.libPaths("lib")
library(metrumrg)## Loading required package: reshape## Loading required package: lattice## Loading required package: XML## Loading required package: MASS## Loading required package: grid## metrumrg 5.57## enter "?metrumrg" for helpsuppressMessages(library(rstan))
suppressMessages(library(bayesplot))
suppressMessages(library(tidyverse))
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinelibrary(parallel)
library(knitr)
library(kableExtra)##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(qapply)## Loading required package: rlecuyer##
## Attaching package: 'qapply'## The following object is masked from 'package:metrumrg':
##
## qstatlibrary(PKPDmisc)
library(GGally)##
## Attaching package: 'GGally'## The following object is masked from 'package:dplyr':
##
## nasa## The following object is masked from 'package:metrumrg':
##
## wrapknitr::opts_chunk$set(
comment = '.',
fig.height = 5,
fig.width = 9,
eval.after = 'fig.cap',
message = FALSE,
fig.path = file.path(figDir, paste(modelName, "_", sep = ""))
)
## Go back to default ggplot2 theme that was overridden by bayesplot
theme_set(theme_gray())
source(file.path(scriptDir, "tools", "stanTools2.R"))
source(file.path(scriptDir, "tools", "functions.R"))
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
RNGkind("L'Ecuyer-CMRG")
set.seed(11191951)
set.seed(10271998) The NONMEM model stub
modelFile <- file.path(modelDir, paste(modelName, "stub.ctl", sep = ""))
modelText <- readLines(modelFile)
cat(modelText, sep = "\n"). $SIZES LVR=50
. $PROBLEM RUN 1 neutropenia PKPD example with Friberg-Karlsson PK-PD
.
. $INPUT C ID TIME EVID AMT CMT DVO DV WT DVID
. $DATA ../../../data/friberg.csv IGNORE=(C='C')
.
.
. $ABBR DECLARE INTEGER FIRST_WRITE_PAR
. $ABBR DECLARE INTEGER FIRST_WRITE_IPAR
.
. $SUBROUTINES ADVAN13 TOL=6 OTHER = ../../extrasend.f90
.
. $MODEL
. NCOMPARTMENTS=8
. COMP = (GUT, DEFDOSE)
. COMP = (CENT)
. COMP = (PERIPH)
. COMP = (PROL)
. COMP = (TRANS1)
. COMP = (TRANS2)
. COMP = (TRANS3)
. COMP = (CIRC)
.
. $PK
. ; Request extra information for Bayesian analysis.
. ; An extra call will then be made for accepted samples
. include '/opt/NONMEM/nm74/util/nonmem_reserved_general'
. BAYES_EXTRA_REQUEST=1
.
. nThin = THETA(12)
.
. VWT = LOG(WT/70) ; normalized to 70 kg adult
.
. MU_1 = THETA(1) + 0.75*VWT ; CL
. MU_2 = THETA(2) + 0.75*VWT ; Q
. MU_3 = THETA(3) + VWT ; V1
. MU_4 = THETA(4) + VWT ; V2
. MU_5 = THETA(5) ; KA
. MU_6 = THETA(6) ; MTT
. MU_7 = THETA(7) ; CIRC0
. MU_8 = THETA(8) ; ALPHA
. MU_9 = THETA(9) ; GAMMA
. MU_10 = THETA(10) ; log(SIGMA) additive error PK
. MU_11 = THETA(11) ; log(SIGMA) additive error PD
.
. " CALL EXTRASEND()
.
. CL = EXP(MU_1 +ETA(1))
. Q = EXP(MU_2 +ETA(2))
. V1 = EXP(MU_3 +ETA(3))
. V2 = EXP(MU_4 +ETA(4))
. KA = EXP(MU_5 +ETA(5))
. MTT = EXP(MU_6 +ETA(6))
. CIRC0=EXP(MU_7 +ETA(7))
. ALPHA=EXP(MU_8 +ETA(8))
.
. GAMMA= EXP(MU_9 + ETA(9))
. SIG1 = EXP(MU_10 + ETA(10)) ; Additive Error PK
. SIG2 = EXP(MU_11 + ETA(11)) ; Additive Error PD
.
. ; intermediate calculations
. K10 = CL/V1
. K12 = Q/V1
. K21 = Q/V2
. KTR = 4/MTT
.
. KPROL = KTR
. KCIRC = KTR
.
. ; initial PK conditions
. A_0(1) = 0
. A_0(2) = 0
. A_0(3) = 0
.
. ; initial PD conditions
. A_0(4) = CIRC0
. A_0(5) = CIRC0
. A_0(6) = CIRC0
. A_0(7) = CIRC0
. A_0(8) = CIRC0
.
.
. S2 = V1
. S3 = V2
.
.
. IF(BAYES_EXTRA==1 .AND. NIREC==1 .AND. NDREC==1 .AND. &
. ITER_REPORT>0 .AND. &
. ITER_REPORT/nThin == INT(ITER_REPORT/nThin)) THEN
. IF(FIRST_WRITE_PAR==0) THEN
. " OPEN(unit=52,FILE='./par.txt')
. FIRST_WRITE_PAR=1
. ENDIF
. " WRITE(52,'(I12,1X,50(1X,1PG19.10E3))') ITER_REPORT, &
. " THETA(1), THETA(2), THETA(3), THETA(4), THETA(5), THETA(6), THETA(7), &
. " THETA(8), THETA(9), THETA(10), THETA(11), &
. " OMEGA(1,1), OMEGA(2,1), OMEGA(2,2), &
. " OMEGA(3,1), OMEGA(3,2), OMEGA(3,3), &
. " OMEGA(4,1), OMEGA(4,2), OMEGA(4,3), OMEGA(4,4), &
. " OMEGA(5,1), OMEGA(5,2), OMEGA(5,3), OMEGA(5,4), OMEGA(5,5), &
. " OMEGA(6,1), OMEGA(6,2), OMEGA(6,3), OMEGA(6,4), OMEGA(6,5), OMEGA(6,6), &
. " OMEGA(7,1), OMEGA(7,2), OMEGA(7,3), OMEGA(7,4), OMEGA(7,5), OMEGA(7,6), OMEGA(7,7), &
. " OMEGA(8,1), OMEGA(8,2), OMEGA(8,3), OMEGA(8,4), OMEGA(8,5), OMEGA(8,6), OMEGA(8,7), OMEGA(8,8)
. ENDIF
.
. IF(BAYES_EXTRA==1 .AND. NDREC==1 .AND. ITER_REPORT>0 .AND. &
. ITER_REPORT/nThin == INT(ITER_REPORT/nThin)) THEN
. IF(FIRST_WRITE_IPAR==0) THEN
. " OPEN(unit=50,FILE='./ipar'//TFI(PNM_NODE_NUMBER)//'.txt')
. FIRST_WRITE_IPAR=1
. ENDIF
. " WRITE(50,'(I12,1X,F14.0,13(1X,1PG12.5))') ITER_REPORT, ID, &
. " ETA(1), ETA(2), ETA(3), ETA(4), ETA(5), ETA(6), ETA(7), ETA(8)
. ENDIF
.
.
.
. $DES
.
. DADT(1) = -KA*A(1)
. DADT(2) = KA*A(1) - (K10 + K12)*A(2) + K21*A(3)
. DADT(3) = K12 *A(2) - K21*A(3)
.
. CONC = A(2)/V1 ; Drug concentration
.
. EDRUG = ALPHA * CONC
. IF(EDRUG > 1.0) EDRUG = 1.0
.
. DADT(4) = KPROL *A(4) *(1 - EDRUG) * ((CIRC0 / A(8)) **GAMMA) - KTR *A(4)
. DADT(5) = KTR * (A(4) - A(5))
. DADT(6) = KTR * (A(5) - A(6))
. DADT(7) = KTR * (A(6) - A(7))
. DADT(8) = KTR*A(7) - KCIRC*A(8)
.
.
. $ERROR (OBSERVATIONS ONLY)
. IPRED = LOG(F)
. IND = 0
. IF(DVID.EQ.2) IND = 1
. YPK = IPRED+SIG1*EPS(1) ; Error for PK
. YPD = IPRED+SIG2*EPS(2) ; Error for PD
. Y = YPK*IND + YPD*(1-IND)
.
. ; Code below into chainsAssemble initial estimates, priors and table specs
## create initial estimates
geninit <- function(){
paste(c(
"; Initial values of THETA",
"$THETA",
paste(rnorm(1, log(10),0.5), " ; Log CL = THETA(1)"),
paste(rnorm(1, log(15), 0.5), " ; Log Q = THETA(2)"),
paste(rnorm(1, log(35), 0.5), " ; Log V1 = THETA(3)"),
paste(rnorm(1, log(105), 0.5), " ; Log V2= THETA(4) "),
paste(rnorm(1, log(2), 0.5), " ; Log Ka = THETA(5) "),
paste(rnorm(1, log(125), 0.2), " ; Log MTT = THETA(6) "),
paste(rnorm(1, log(5), 0.2), " ; Log CIRC0 = THETA(7) "),
paste(rnorm(1, log(3E-4), 1), " ; Log ALPHA = THETA(8) "),
paste(rnorm(1, log(0.17), 0.2), " ; Log GAMMA = THETA(9) "),
paste(rnorm(1, log(0.2), 0.2), " ; Log PK Add Error = THETA(10) "),
paste(rnorm(1, log(0.2), 0.2), " ; Log PD Add Error = THETA(11) "),
"$OMEGA BLOCK(8) ;INITIAL values of OMEGAs",
paste(exp(rnorm(1, log(0.25), 0.5))),
paste(rep(rnorm(1, 0.001, 0.01),1)), paste(exp(rnorm(1, log(0.25), 0.5))),
paste(rep(rnorm(1, 0.001, 0.01),2)), paste(exp(rnorm(1, log(0.25), 0.5))),
paste(rep(rnorm(1, 0.001, 0.01),3)), paste(exp(rnorm(1, log(0.25), 0.5))),
paste(rep(rnorm(1, 0.001, 0.01),4)), paste(exp(rnorm(1, log(0.25), 0.5))),
paste(rep(rnorm(1, 0.001, 0.01),5)), paste(exp(rnorm(1, log(0.25), 0.5))),
paste(rep(rnorm(1, 0.001, 0.01),6)), paste(exp(rnorm(1, log(0.25), 0.5))),
paste(rep(rnorm(1, 0.001, 0.01),7)), paste(exp(rnorm(1, log(0.05), 0.5))),
"$OMEGA DIAG(3) ;INITIAL values of unused OMEGAs",
"0 FIX; GAMMA",
"0 FIX; Add PK",
"0 FIX; Add PD",
"$SIGMA ;Initial value of SIGMA",
"1 FIX; Add PK",
"1 FIX; Add PD"),
sep = "\n")
}
## Set parameters of the prior distribution
priors <- paste(c("$PRIOR NWPRI",
"$THETAP ; Prior information of THETAS",
"(2.3 FIX) ; CL Log(10) ",
"(2.7 FIX) ; Q Log(15)" ,
"(3.56 FIX) ; V1 Log(35)",
"(4.7 FIX) ; V2 Log(105)",
"(0.693 FIX) ; KA log(2)",
"(4.83 FIX) ; MTT Log(125)",
"(1.61 FIX) ; CIRC0 Log(5) ",
"(-8.11 FIX) ; ALPHA Log(3E-4)",
"(-1.78 FIX) ; GAMMA Log(0.17)",
"(-1.60 FIX) ; log Add error PK Log(0.20)",
"(-1.60 FIX) ; log Add error PD Log(0.20)",
"$THETAPV BLOCK(11) ; variances for priors on THETAS (var-cov)",
"0.25 FIX ; CL ",
"0.00 0.25 ; Q",
"0.00 0.00 0.25 ; V1",
"0.00 0.00 0.00 0.25 ; V2",
"0.00 0.00 0.00 0.00 0.25 ; KA ",
"0.00 0.00 0.00 0.00 0.00 0.04 ; MTT ",
"0.00 0.00 0.00 0.00 0.00 0.00 0.04 ; CIRC0 ",
"0.00 0.00 0.00 0.00 0.00 0.00 0.00 1 ; ALPHA ",
"0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.04 ; GAMMA ",
"0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 ; Error PK ",
"0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 ; Error PD ",
"$OMEGAP BLOCK(8) ;prior information for OMEGA",
"0.045 FIX ; CL ",
"0.00 0.045 ; Q",
"0.00 0.00 0.045 ; V1",
"0.00 0.00 0.00 0.045 ; V2",
"0.00 0.00 0.00 0.00 0.045 ; KA",
"0.00 0.00 0.00 0.00 0.00 0.045 ; MTT",
"0.00 0.00 0.00 0.00 0.00 0.00 0.045 ; CIRC0",
"0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.045 ; ALPHA",
"$OMEGAPD (11, FIXED) ; df for OMEGA prior"
),
sep = "\n")
## Table statements
tables <- paste(c("$TABLE ID EVID TIME DV IPRED CWRES CWRESI NPDE WT",
"NOPRINT ONEHEADER FILE=./.tab",
"$TABLE ID WT CL Q V1 V2 KA MTT CIRC0 ALPHA ETAS(1:LAST)",
"NOPRINT ONEHEADER FILE=./par.tab"),
sep = "\n")Run NONMEM
nChains <- 4
nPost <- 250 ## Number of post-burn-in samples per chain after thinning
nBurn <- 250 ## Number of burn-in samples per chain after thinning
nThin <- 1
seed = sample(10000, size = nChains)
if(fitModel){
if(!file.exists(file.path(modelDir, modelName)))
dir.create(file.path(modelDir, modelName))
runChains <- lapply(1:nChains, runChain,
modelName = modelName, modelDir = modelDir,
priors = priors, tables = tables,
nPost = nPost, nBurn = nBurn, nThin = nThin,
seed = seed, print = 100,
OLKJDF = 0, OVARF = 1, AUTO = 2,
NUTS_DELTA = 0.8,
grid = TRUE, method = "NUTS", mode='not para')
}read in data file
dataFile <- "friberg.csv"
nmData <- read.csv(file.path(dataDir, dataFile), as.is = TRUE, na.strings = '.')Get population parameters
chains <- 1:nChains
##chains <- c(1, 3, 4)
nChains2 <- length(chains)
## Read in posterior distributions by chain and add a column called chain
popParameters = map(chains, function(thisChain){
param <- data.table::fread(file.path(modelDir, modelName, paste(modelName, ".", thisChain, sep = ""),
"par.txt"), data.table = FALSE)
names(param) <- c("iteration",
"THETA1", "THETA2", "THETA3", "THETA4",
"THETA5", "THETA6", "THETA7", "THETA8",
"THETA9", "THETA10", "THETA11",
"OMEGA11", "OMEGA21", "OMEGA22","OMEGA13", "OMEGA23", "OMEGA33",
"OMEGA14", "OMEGA24", "OMEGA34","OMEGA44", "OMEGA15", "OMEGA25",
"OMEGA35", "OMEGA45", "OMEGA55","OMEGA16", "OMEGA26", "OMEGA36",
"OMEGA46", "OMEGA56", "OMEGA66","OMEGA17", "OMEGA27", "OMEGA37",
"OMEGA47", "OMEGA57", "OMEGA67","OMEGA77", "OMEGA18", "OMEGA28",
"OMEGA38", "OMEGA48", "OMEGA58","OMEGA68", "OMEGA78", "OMEGA88")
param %>% mutate(chain = thisChain)
}) %>% bind_rows %>%
mutate(sample = 1:n())
## Calculate more interpretable paraneters
popParameters2 <- popParameters %>%
mutate(CLHat = exp(THETA1),
QHat = exp(THETA2),
V1Hat = exp(THETA3),
V2Hat = exp(THETA4),
KaHat = exp(THETA5),
MTTHat = exp(THETA6),
CIRC0Hat = exp(THETA7),
ALPHAHat = exp(THETA8),
GAMMAHat=exp(THETA9),
SigmaPK = exp(THETA10),
SigmaPD = exp(THETA11)) %>%
select(starts_with("OMEGA"),
CLHat, QHat, V1Hat, V2Hat, KaHat,
MTTHat, CIRC0Hat, ALPHAHat, GAMMAHat, SigmaPK,
SigmaPD, chain, iteration, sample)
## Convert to 3-D array with dims = {iterations, chains, parameters}.
popParArray <- array(double(nPost * nChains2 * (ncol(popParameters2) - 3)),
dim = c(nPost, nChains2, ncol(popParameters2) - 3),
dimnames = list(NULL, NULL, setdiff(names(popParameters2), c("chain", "iteration", "sample"))))
for(iChain in 1:nChains2){
popParArray[,iChain,] <- popParameters2 %>%
filter(chain == chains[iChain],
iteration > 0) %>%
select(-chain, -iteration, -sample) %>%
as.matrix
}Get individual parameters
indParameters = map(chains, function(thisChain){
files <- list.files(file.path(modelDir, modelName, paste(modelName, ".", thisChain, sep = "")))
iparFiles <- files[grepl("^ipar", files)]
param <- map(iparFiles, function(thisFile){
data.table::fread(file.path(modelDir, modelName, paste(modelName, ".", thisChain, sep = ""),
thisFile), data.table = FALSE)
}) %>%
bind_rows %>% select(1:10)
names(param) <- c("iteration", "ID", "ETA1", "ETA2",
"ETA3","ETA4","ETA5","ETA6","ETA7","ETA8")
param %>% mutate(chain = thisChain)
}) %>% bind_rows %>%
arrange(chain, iteration, ID) %>%
left_join(popParameters)MCMC diagnostics and posterior distributions of parameters
options(bayesplot.base_size = 12,
bayesplot.base_family = "sans")
color_scheme_set(scheme = "brightblue")
myTheme <- theme(text = element_text(size = 12),
axis.text = element_text(size = 12))
plot_mcmcHistory <- mcmcHistory(popParArray,
nParPerPage = 5, myTheme = myTheme)
plot_mcmcDensityByChain <- mcmcDensity(popParArray,
nParPerPage = 16, byChain = TRUE,
myTheme = theme(text = element_text(size = 12),
axis.text = element_text(size = 10)))
plot_mcmcDensity <- mcmcDensity(popParArray, nParPerPage = 16,
myTheme = theme(text = element_text(size = 12),
axis.text = element_text(size = 10)))
plot_mcmcHistory. [[1]]
Friberg-Karlsson PK-PD example: MCMC trace plots.
.
. [[2]]
Friberg-Karlsson PK-PD example: MCMC trace plots.
.
. [[3]]
Friberg-Karlsson PK-PD example: MCMC trace plots.
.
. [[4]]
Friberg-Karlsson PK-PD example: MCMC trace plots.
.
. [[5]]
Friberg-Karlsson PK-PD example: MCMC trace plots.
.
. [[6]]
Friberg-Karlsson PK-PD example: MCMC trace plots.
.
. [[7]]
Friberg-Karlsson PK-PD example: MCMC trace plots.
.
. [[8]]
Friberg-Karlsson PK-PD example: MCMC trace plots.
.
. [[9]]
Friberg-Karlsson PK-PD example: MCMC trace plots.
.
. [[10]]
Friberg-Karlsson PK-PD example: MCMC trace plots.
plot_mcmcDensityByChain. [[1]]
Friberg-Karlsson PK-PD example: Univariate marginal densities by chain.
.
. [[2]]
Friberg-Karlsson PK-PD example: Univariate marginal densities by chain.
.
. [[3]]
Friberg-Karlsson PK-PD example: Univariate marginal densities by chain.
plot_mcmcDensity. [[1]]
Friberg-Karlsson PK-PD example: Univariate marginal densities.
.
. [[2]]
Friberg-Karlsson PK-PD example: Univariate marginal densities.
.
. [[3]]
Friberg-Karlsson PK-PD example: Univariate marginal densities.
#Not ideal when there are many parameters
#mcmc_pairs(popParArray[,,!grepl("^OMEGA", dimnames(popParArray)[[3]])])
#mcmc_pairs(popParArray[,,grepl("^OMEGA", dimnames(popParArray)[[3]])])
ggpairs(filter(popParameters2) %>%
select(CLHat:GAMMAHat),
diag = list(continuous = 'barDiag'),
upper = list(continuous = 'points'),
lower=list(continuous='points'))
Friberg-Karlsson PK-PD example: Correlation Plots.
ggpairs(filter(popParameters) %>%
select(OMEGA11,OMEGA22,OMEGA33,OMEGA44,
OMEGA55,OMEGA66,OMEGA77,OMEGA88),
diag = list(continuous = 'barDiag'),
upper = list(continuous = 'points'),
lower=list(continuous='points'))
Friberg-Karlsson PK-PD example: Correlation Plots.
caption <- paste("Friberg-Karlsson PK-PD example:",
c(rep("MCMC trace plots.", length(plot_mcmcHistory)),
rep("Univariate marginal densities by chain.",
length(plot_mcmcDensityByChain)),
rep("Univariate marginal densities.",
length(plot_mcmcDensityByChain)),
rep("Correlation Plots.", 2)))
ptable <- monitor(popParArray,
warmup = 0, print = FALSE) %>%
as.matrix() %>%
formatC(3) %>%
as.data.frame
write.csv(ptable, file = file.path(tabDir, paste(modelName,
"ParameterTable.csv", sep = "")))
ptable %>%
rename(SEmean = se_mean, SD = sd, pct2.5 = "2.5%", pct25 = "25%", median = "50%",
pct75 = "75%", pct97.5 = "97.5%", Neff = "n_eff") %>%
mutate(parameter = rownames(.), "95% CI" = paste("(", pct2.5, ", ", pct97.5, ")",
sep = "")) %>%
select(parameter, mean, SD, median, "95% CI", Neff, Rhat) %>%
kable(caption = "Summary of model parameter estimates.") %>%
kable_styling(bootstrap_options = "striped", full_width = F)| parameter | mean | SD | median | 95% CI | Neff | Rhat |
|---|---|---|---|---|---|---|
| OMEGA11 | 0.135 | 0.0866 | 0.111 | (0.0434, 0.371) | 652 | 1 |
| OMEGA21 | -0.0336 | 0.0667 | -0.027 | (-0.191, 0.0861) | 840 | 0.999 |
| OMEGA22 | 0.169 | 0.104 | 0.142 | (0.0546, 0.479) | 597 | 1 |
| OMEGA13 | -0.00534 | 0.062 | -0.00315 | (-0.133, 0.106) | 553 | 1 |
| OMEGA23 | -0.00751 | 0.0669 | -0.00681 | (-0.151, 0.13) | 621 | 1 |
| OMEGA33 | 0.11 | 0.0761 | 0.0866 | (0.0359, 0.329) | 713 | 1 |
| OMEGA14 | 0.0525 | 0.0782 | 0.0384 | (-0.0632, 0.25) | 692 | 1 |
| OMEGA24 | -0.0724 | 0.09 | -0.0541 | (-0.31, 0.0596) | 673 | 1 |
| OMEGA34 | -0.00232 | 0.0749 | 0.000888 | (-0.168, 0.141) | 725 | 0.999 |
| OMEGA44 | 0.199 | 0.142 | 0.162 | (0.0611, 0.585) | 715 | 1 |
| OMEGA15 | -0.00232 | 0.0662 | -0.00181 | (-0.13, 0.12) | 614 | 1 |
| OMEGA25 | 0.00947 | 0.0646 | 0.00845 | (-0.128, 0.153) | 451 | 1 |
| OMEGA35 | 0.0196 | 0.0603 | 0.0121 | (-0.0775, 0.158) | 753 | 1 |
| OMEGA45 | -0.0165 | 0.0776 | -0.0132 | (-0.18, 0.122) | 532 | 1 |
| OMEGA55 | 0.117 | 0.0833 | 0.0938 | (0.0353, 0.32) | 878 | 1 |
| OMEGA16 | -0.016 | 0.0526 | -0.0139 | (-0.12, 0.0849) | 578 | 1 |
| OMEGA26 | 0.00682 | 0.0558 | 0.00827 | (-0.116, 0.124) | 613 | 1 |
| OMEGA36 | 0.00529 | 0.047 | 0.00345 | (-0.0756, 0.105) | 737 | 1 |
| OMEGA46 | -0.00318 | 0.0619 | -0.00595 | (-0.118, 0.135) | 698 | 1.01 |
| OMEGA56 | -0.00348 | 0.0509 | -0.00131 | (-0.109, 0.0931) | 891 | 1 |
| OMEGA66 | 0.095 | 0.0668 | 0.0771 | (0.0339, 0.277) | 673 | 1 |
| OMEGA17 | 0.0144 | 0.0574 | 0.0126 | (-0.0901, 0.143) | 684 | 1 |
| OMEGA27 | 0.0354 | 0.0666 | 0.0255 | (-0.0659, 0.199) | 603 | 1.01 |
| OMEGA37 | -0.00556 | 0.0537 | -0.00447 | (-0.11, 0.101) | 706 | 1 |
| OMEGA47 | -0.0211 | 0.0764 | -0.0157 | (-0.181, 0.1) | 773 | 1.01 |
| OMEGA57 | 0.0104 | 0.0552 | 0.00674 | (-0.0949, 0.131) | 624 | 1 |
| OMEGA67 | -0.014 | 0.0483 | -0.011 | (-0.108, 0.0711) | 713 | 1 |
| OMEGA77 | 0.121 | 0.0815 | 0.101 | (0.0412, 0.339) | 493 | 1.01 |
| OMEGA18 | 0.0372 | 0.062 | 0.0287 | (-0.0583, 0.167) | 586 | 1 |
| OMEGA28 | -0.0206 | 0.0628 | -0.0165 | (-0.159, 0.104) | 777 | 1 |
| OMEGA38 | 0.00601 | 0.052 | 0.0029 | (-0.0856, 0.124) | 669 | 1 |
| OMEGA48 | 0.0265 | 0.0672 | 0.0212 | (-0.0977, 0.17) | 662 | 1 |
| OMEGA58 | -0.00619 | 0.0583 | -0.00455 | (-0.132, 0.114) | 759 | 1.01 |
| OMEGA68 | -0.0103 | 0.0488 | -0.00807 | (-0.0983, 0.0716) | 427 | 1 |
| OMEGA78 | 0.0148 | 0.0564 | 0.0117 | (-0.0807, 0.134) | 497 | 1 |
| OMEGA88 | 0.122 | 0.0836 | 0.101 | (0.0421, 0.34) | 531 | 1.01 |
| CLHat | 11.2 | 1.57 | 11 | (8.35, 14.6) | 832 | 1 |
| QHat | 22.6 | 3.63 | 22.4 | (15.9, 30.2) | 792 | 1 |
| V1Hat | 29.2 | 5.51 | 29.1 | (18.7, 40.5) | 801 | 1 |
| V2Hat | 117 | 20.4 | 114 | (82.9, 162) | 833 | 0.999 |
| KaHat | 1.69 | 0.344 | 1.66 | ( 1.1, 2.42) | 712 | 1 |
| MTTHat | 118 | 13.1 | 118 | (93.7, 145) | 741 | 1 |
| CIRC0Hat | 5.75 | 0.686 | 5.73 | (4.46, 7.15) | 956 | 1 |
| ALPHAHat | 0.000323 | 4.72e-05 | 0.000321 | (0.000236, 0.000426) | 780 | 1 |
| GAMMAHat | 0.151 | 0.013 | 0.15 | (0.127, 0.177) | 1.05e+03 | 1.01 |
| SigmaPK | 0.1 | 0.00442 | 0.1 | (0.092, 0.109) | 1.06e+03 | 1 |
| SigmaPD | 0.106 | 0.00678 | 0.106 | (0.0941, 0.121) | 908 | 1 |
Construct design data frame
## For now just use the NONMEM data set
design <- nmData %>%
mutate(ID = id, WT = weight) %>%
filter(is.na(C))mrgsolve model code for posterior predictions
library(mrgsolve)
indCode =
'
$PARAM
THETA1 = 2.3
THETA2 = 2.7
THETA3 = 3.56
THETA4 = 4.7
THETA5 = 0.693
THETA6 = 4.43
THETA7 = 2.33
THETA8 = -1.78
THETA9 = -8.11
THETA10 = -1
THETA11 = -1
WT = 70
ETA1 = 0
ETA2 = 0
ETA3 = 0
ETA4 = 0
ETA5 = 0
ETA6 = 0
ETA7 = 0
ETA8 = 0
$CMT GUT, CENT, PERIPH,
PROL, TRANS1, TRANS2, TRANS3, CIRC
$MAIN
double VWT = log(WT/70);
double CL = exp(THETA1 + (VWT * 0.75) + ETA1 + ETA(1));
double Q = exp(THETA2 + (VWT * 0.75) + ETA2 + ETA(2));
double V1 = exp(THETA3 + VWT + ETA3 + ETA(3));
double V2 = exp(THETA4 + VWT + ETA4 + ETA(4));
double KA = exp(THETA5 + ETA5 + ETA(5));
double MTT = exp(THETA6 + ETA6 + ETA(6));
double CIRC0 = exp(THETA7 + ETA7 + ETA(7));
double ALPHA = exp(THETA8 + ETA8 + ETA(8));
double GAMMA = exp(THETA9);
double SIG1 = exp(THETA10);
double SIG2 = exp(THETA11);
// intermediate calculations
double K10 = CL/V1;
double K12 = Q/V1;
double K21 = Q/V2;
double KTR = 4/MTT;
double KPROL = KTR;
double KCIRC = KTR;
PROL_0 = CIRC0;
TRANS1_0 = CIRC0;
TRANS2_0 = CIRC0;
TRANS3_0 = CIRC0;
CIRC_0 = CIRC0;
$ODE
dxdt_GUT = -KA*GUT;
dxdt_CENT = KA*GUT - K10*CENT - K12*CENT + K21*PERIPH;
dxdt_PERIPH = K12*CENT - K21*PERIPH;
double CONC = std::max(0.0,CENT/V1);
double EDRUG = CONC*ALPHA;
if(EDRUG > 1) EDRUG = 1;
dxdt_PROL = KPROL* PROL *(1 - EDRUG) * (pow((CIRC0 / CIRC), GAMMA)) - KTR *PROL;
dxdt_TRANS1 = KTR * (PROL - TRANS1);
dxdt_TRANS2 = KTR * (TRANS1 - TRANS2);
dxdt_TRANS3 = KTR * (TRANS2 - TRANS3);
dxdt_CIRC = KTR*TRANS3 - KCIRC*CIRC;
$SIGMA 1 1
$OMEGA @block
0
0 0
0 0 0
0 0 0 0
0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
$TABLE
double cPredHat = std::max(0.000001,CENT/V1);
double lcPred = log(cPredHat) + SIG1*EPS(1);
double cPred = exp(lcPred);
double nPredHat = CIRC;
double lnPred = log(CIRC) + SIG2*EPS(2);
double nPred = exp(lnPred);
$CAPTURE lcPred lnPred cPred nPred
'
## Compile model
friMod <- mcode("FRImodel", indCode)
if(FALSE){
## Test runs
out <- friMod %>%
data_set(design) %>%
carry.out(DV) %>%
mrgsim %>% as.data.frame()
ggplot(data=out, aes(x=time, y=cPred)) + geom_point(aes(group=ID, color=as.factor(ID)))
}Individual predictions
samples <- unique(indParameters$sample)
set.seed(10301987)
indSim1 <- mclapply(samples,
function(thisSample){
friMod %>%
data_set(design) %>%
idata_set(indParameters %>%
filter(sample == thisSample)) %>%
carry_out(sample, ID, time, DV, evid, DVID) %>%
mrgsim() %>%
as.tbl
}) %>%
bind_rowsPopulation predictions
samples <- unique(popParameters$sample)
set.seed(10301987)
popSim1 <- mclapply(samples,
function(thisSample){
friMod %>%
data_set(design) %>%
param(popParameters %>%
filter(sample == thisSample)) %>%
omat(omat(as_bmat(popParameters %>%
filter(sample == thisSample),"OMEGA"))) %>%
carry_out(sample, ID, time, DV, evid, DVID) %>%
mrgsim %>%
as.tbl
}) %>%
bind_rowsPosterior predictive distributions
## Combined plots
indPredC <- indSim1 %>% filter(evid==0, DVID==2) %>%
group_by(ID, time, evid) %>%
summarize(lbInd = quantile(cPred, probs = 0.05, na.rm = TRUE),
medianInd = quantile(cPred, probs = 0.5, na.rm = TRUE),
ubInd = quantile(cPred, probs = 0.95, na.rm = TRUE))
indPredN <- indSim1 %>% filter(evid==0, DVID==1) %>%
group_by(ID, time, evid) %>%
summarize(lbInd = quantile(nPred, probs = 0.05, na.rm = TRUE),
medianInd = quantile(nPred, probs = 0.5, na.rm = TRUE),
ubInd = quantile(nPred, probs = 0.95, na.rm = TRUE))
popPredC <- popSim1 %>% filter(evid==0, DVID==2) %>%
group_by(ID, time, evid) %>%
summarize(lbPop = quantile(cPred, probs = 0.05, na.rm = TRUE),
medianPop = quantile(cPred, probs = 0.5, na.rm = TRUE),
ubPop = quantile(cPred, probs = 0.95, na.rm = TRUE))
popPredN <- popSim1 %>% filter(evid==0, DVID==1) %>%
group_by(ID, time, evid) %>%
summarize(lbPop = quantile(nPred, probs = 0.05, na.rm = TRUE),
medianPop = quantile(nPred, probs = 0.5, na.rm = TRUE),
ubPop = quantile(nPred, probs = 0.95, na.rm = TRUE))
predAllC <- indPredC %>%
left_join(popPredC) %>%
left_join(design %>% filter(evid == 0, DVID==2)) %>%
mutate(DV = suppressWarnings(as.numeric(DV)))
predAllN <- indPredN %>%
left_join(popPredN) %>%
left_join(design %>% filter(evid == 0, DVID==1)) %>%
mutate(DV = suppressWarnings(as.numeric(DV)))
p1 <- ggplot(predAllC, aes(x = time, y = DV))
p2 <- p1 +
geom_line(aes(x = time, y = medianPop,
color = "population")) +
geom_ribbon(aes(ymin = lbPop, ymax = ubPop,
fill = "population"),
alpha = 0.25) +
geom_line(aes(x = time, y = medianInd,
color = "individual")) +
geom_ribbon(aes(ymin = lbInd, ymax = ubInd,
fill = "individual"),
alpha = 0.25) +
scale_color_brewer(name ="prediction",
breaks=c("individual", "population"),
palette = "Set1") +
scale_fill_brewer(name ="prediction",
breaks=c("individual", "population"),
palette = "Set1")
p3 <- p2 + geom_point() +
labs(title = 'PK Model',
x = "Time (h)",
y = "PK Conc (ng/mL)") +
theme(text = element_text(size = 12),
axis.text = element_text(size = 12),
## legend.position = c(0.8, 0.25),
strip.text = element_text(size = 8)) +
facet_wrap(~ ID, scales = "free_y")
####PD Plots#####
p1b <- ggplot(predAllN, aes(x = time, y = DV))
p2b <- p1b +
geom_line(aes(x = time, y = medianPop,
color = "population")) +
geom_ribbon(aes(ymin = lbPop, ymax = ubPop,
fill = "population"),
alpha = 0.25) +
geom_line(aes(x = time, y = medianInd,
color = "individual")) +
geom_ribbon(aes(ymin = lbInd, ymax = ubInd,
fill = "individual"),
alpha = 0.25) +
scale_color_brewer(name ="prediction",
breaks=c("individual", "population"),
palette = "Set1") +
scale_fill_brewer(name ="prediction",
breaks=c("individual", "population"),
palette = "Set1")
p3b <- p2b + geom_point() +
labs(title = 'PD Model',
x = "Time (h)",
y = "ANC (10^9/L)") +
theme(text = element_text(size = 12),
axis.text = element_text(size = 12),
## legend.position = c(0.8, 0.25),
strip.text = element_text(size = 8)) +
facet_wrap(~ ID, scales = "free_y")
####
p3
p3b
sessionInfo(). R version 3.5.1 (2018-07-02)
. Platform: x86_64-pc-linux-gnu (64-bit)
. Running under: Ubuntu 14.04.5 LTS
.
. Matrix products: default
. BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
. LAPACK: /usr/lib/lapack/liblapack.so.3.0
.
. locale:
. [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
. [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
. [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
. [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
. [9] LC_ADDRESS=C LC_TELEPHONE=C
. [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
.
. attached base packages:
. [1] parallel grid stats graphics grDevices utils datasets
. [8] methods base
.
. other attached packages:
. [1] mrgsolve_0.9.2 GGally_1.4.0 PKPDmisc_2.1.1
. [4] qapply_1.39 rlecuyer_0.3-4 kableExtra_1.1.0
. [7] knitr_1.23 gridExtra_2.3 forcats_0.4.0
. [10] stringr_1.4.0 dplyr_0.8.3 purrr_0.3.2
. [13] readr_1.3.1 tidyr_0.8.3 tibble_2.1.3
. [16] tidyverse_1.2.1 bayesplot_1.7.0 rstan_2.19.2
. [19] ggplot2_3.2.0 StanHeaders_2.18.1-10 metrumrg_5.57
. [22] MASS_7.3-51.1 XML_3.98-1.20 lattice_0.20-38
. [25] reshape_0.8.8
.
. loaded via a namespace (and not attached):
. [1] httr_1.4.0 jsonlite_1.6
. [3] viridisLite_0.3.0 modelr_0.1.4
. [5] assertthat_0.2.1 highr_0.8
. [7] stats4_3.5.1 cellranger_1.1.0
. [9] yaml_2.2.0 pillar_1.4.2
. [11] backports_1.1.4 glue_1.3.1
. [13] digest_0.6.20 RColorBrewer_1.1-2
. [15] rvest_0.3.4 colorspace_1.4-1
. [17] htmltools_0.3.6 plyr_1.8.4
. [19] pkgconfig_2.0.2 broom_0.5.2
. [21] haven_2.1.1 scales_1.0.0
. [23] webshot_0.5.1 processx_3.4.1
. [25] generics_0.0.2 withr_2.1.2
. [27] lazyeval_0.2.2 cli_1.1.0
. [29] magrittr_1.5 crayon_1.3.4
. [31] readxl_1.3.1 evaluate_0.14
. [33] ps_1.3.0 nlme_3.1-137
. [35] RcppArmadillo_0.9.600.4.0 xml2_1.2.1
. [37] pkgbuild_1.0.3 data.table_1.12.2
. [39] tools_3.5.1 loo_2.1.0
. [41] prettyunits_1.0.2 hms_0.5.0
. [43] matrixStats_0.54.0 munsell_0.5.0
. [45] callr_3.3.1 compiler_3.5.1
. [47] rlang_0.4.0 ggridges_0.5.1
. [49] rstudioapi_0.10 labeling_0.3
. [51] rmarkdown_1.14 gtable_0.3.0
. [53] inline_0.3.15 reshape2_1.4.3
. [55] R6_2.4.0 lubridate_1.7.4
. [57] zeallot_0.1.0 stringi_1.4.3
. [59] Rcpp_1.0.2 vctrs_0.2.0
. [61] tidyselect_0.2.5 xfun_0.8